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Abstract 

We report results of systematic numerical studies of 2D matter-wave soliton families supported 
by an external potential, in a vicinity of the junction between stable and unstable branches of the 
families, where the norm of the solution attains a minimum, facilitating the creation of the soliton. 
The model is based on the Gross-Pitaevskii equation for the self-attractive condensate loaded into 
a quasiperiodic (QP) optical lattice (OL). The same model applies to spatial optical solitons in QP 
photonic crystals. Dynamical properties and stability of the solitons are analyzed with respect to 
variations of the depth and wavenumber of the OL. In particular, it is found that the single-peak 
solitons are stable or not in exact accordance with the Vakhitov-Kolokolov (VK) criterion, while 
double-peak solitons, which are found if the OL wavenumber is small enough, are always unstable 
against splitting. 
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Introduction and the model. A challenging subject in studies of dynamical patterns in 
Bose-Einstein condensates (BECs) and nonlinear optics is the creation of matter-wave or 
photonic solitons in multidimensional settings [l-3|. Various routes to the making of stable 
two- and three-dimensional (2D and 3D) fundamental and vortical solitons have been elab- 
orated theoretically. As demonstrated in Refs. j4|-j^, universal stabilization methods for 
the matter-wave and optical solitons are provided, respectively, by optical lattices (OLs) or 
photonic crystals, i.e., essentially, by spatially periodic potentials. OLs are induced, as inter- 
ference patterns, by coherent laser beams illuminating the condensate in opposite directions, 
while photonic lattices may be created, by means of various technologies, as permanent struc- 
;ures in optical waveguides, or as virtual photoinduced structures in photorefractive crystals 
2! . A more difficult but also realistic possibility is stabilizing solitons by means of nonlinear 
,att.ce, .e.. spatially penodic o,oduMio,. of the nonUneanty coefficient g In principle, 
similar methods may be applied to a gas of polaritons [7|, where the evidence of the BEG 
state was reported too [sl, using properly engineered superlattices [q]. 

The stabilization of 2D and 3D solitons is possible with the help of the fully-dimensional 
OL, whose dimension is equal to that of the entire space, D, and by low- dimensional lattices. 



with dimension D 



1 flQ, 



Other methods 



br the creation of robust solitons rely 
on the t,n,e-penod. ».„...»e„. M of nonHnea. flQ o. Hnea. Q cWcten... of 
the condensate (following the method proposed 17| and later implemented experimentally 
[isl l for the stabilization of 2D solitons in optics by means of the periodic modulation of 
the Kerr coefficient along the propagation distance ). In these contexts, the stability of the 
matter waves in 2D OLs, and under various scenarios of the time-periodic management, 
has been studied extensively, see, e.g., Refs. [l9|, l20|. In addition, the stabilization of 



multidimensional solitons may be provided by nonlocal (dipole-dipole) interaction between 



atoms 



2l| or nonlocal (thermal) nonlinearity in optics |22| . 



Besides periodic OLs, quasiperiodic (QP) ones have also drawn a great deal of interest 
in particular, as the simplest setting for the realization of the Anderson localization o 



matter waves 23|]. The self-trapping of 2D solitons in QP potentials was studied too 
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251]. The objective of this work is to extend the previously reported analysis of the 



stabilization of 2D solitons by lattice potentials to the case of QP lattices and self-attractive 
nonlinearity (negative scattering length of inter-atomic interactions in the BEG), which can 
be readily implemented in ^Li and ®^Rb condensates [26|, and corresponds to the usual 
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Kerr nonlinearity in optics. As known from the previous analyses , UM ; dependence 
between the chemical potential and the norm (which is proportional to the number of atoms 
in BEC, or total power of the optical beam) for 2D solitons supported by lattice potentials, 
IJ.{N), features two branches, stable and unstable ones [with dfi/dN < and dfi/dN > 0, 
respectively, according to the Vakhitov-Kolokolov (VK) criterion [27]]. The branches merge 
at a threshold (minimal) value of A^, below which the solitons decay due to the delocalization 
transition 

Our analysis is based on the 2D Gross-Pitaevskii equation for the BEC mean-field wave 
function, (x, y, t), written in the dimensionless form assuming the self-attractive nonlin- 



earity 



30|: 



where the QP lattice potential of depth 2Vo is taken as jsl, [2^ |25] 



M 



-\/(x,2/) = -K)$^cos(k(")r), (2) 



n=l 



with the set of wave vectors k^") = Hcos(27r(n - 1)/M) , sin (27r(n - 1)/M)} and M = 5 



or M > 7. Here, following Ref. 2^, we focus on the basic case of the Penrose-tiling 
potential, corresponding to M = 5. The 2D profile of the potential is displayed below in 
Fig. [3t^d). Setting Vq > 0, the center of the 2D soliton will be placed at the local minimum 
of potential i^, x = y = 0. The solitons will be characterized by their norm, defined as 
usual: N = j J \'${x,y)f dxdy. The relation of to the actual number of atoms in the 
condensate, Af, is given by means of standard rescaling [30]: Af = {a^/A'nas)N, where a± 
(typically, ~ yum) and (~ 0.1 nm) are the transverse trapping length of the condensate 
and scattering length of the atomic collisions, respectively. In optics, the same equation ([1]), 
with t replaced by the propagation distance, z, governs, the transmission of electromagnetic 
waves with local amplitude ^ in the bulk waveguide with the transverse QP modulation of 
the refractive index. In the latter case, A^ is proportional to the beam's total power. 

Numerical results: soliton families. Simulations of Eq. ([T]) were performed on the 2D 
numerical grid of size 128 x 128, starting with the input in the form of an isotropic Gaussian, 

^'(a;,y) = Aoexp(-g(x' + y2)). (3) 
Initial amplitude Aq, along with the OL depth and wavenumber, Vq and fc, were varied. 
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while the initial width was fixed by setting g = 0.9 [which is possible by means of rescaling 
ofEq. m- 

Before proceeding to numerical results, it is relevant to note that, although the application 
of the variational approximation, which is a ubiquitous analytical tool for the study of bound 
states in nonlinear systems jl, 3|, to 2D solitons in QP potentials is possible {g], the simplest 
isotropic ansatz, taken in the same form as Gaussian ([3]), cannot capture peculiarities of 
the setting based on the QP potential. Indeed, the part of the Lagrangian accounting 
for the interaction of ansatz (j3]) with the underlying OL potential (|2]) consists of integrals 
like VqAIJ J cos (k^^^r) exp(— 2gr^)dr = tt[Vo/ (2g)] exp [— /c^/ (8g)]. Being insensitive to 
the particular orientation of wave vectors k*^"'^ this approximation is too coarse. It may 
be improved by using an anisotropic ansatz, but this will render the variational analysis 
cumbersome. 

The first objective is to construct families of localized ground-state modes, in the form 
of \E'(x,?/,t) = exp{—ifit)ip{x,y), with real wave function {p{x,y) found by means of the 
accelerated imaginary-time method 3l|]- Following the convention commonly adopted in 
physics literature Q-Q, 10|-[l5|, we refer to these modes as "solitons", even though they 
do not feature the unhindered motion characteristic to "genuine" solitons. The simulations 
of Eq. ([1]), rewritten in the imaginary time with a fixed value of /i, quickly converge to the 
ground state, with < 1000 iterations necessary to reduce the residual error to the level of 
10-10. 

In Fig. [H chemical potential /i of the ground state is shown, as a function of its norm 
A^, at two fixed wavenumbers of the Penrose-tiling potential, = 1 (a) and k = 1.5 (b) 
and various values of its depth, Vq. Further, Fig. |2] shows fi{N) for fixed Vq and different 
values of k. Labels Cj and Aj (j = 1,2,3,4) indicate branches which are expected to be 



stable and unstable according to the Vakhitov-Kolokolov (VK) criterion [27|, |28|, i.e., with 
dfi/dN < and dfi/dN > 0, respectively. The imaginary-time algorithm, which generated 
the solitons, ceased to converge at lower termination points of the branches shown in Figs. 
[Hand [21 where the amplitude of the solution becomes too large. 

Points Bj in Figs. [1] and [2] mark the junctions between the stable and unstable branches, 
where dfi/dN diverges, while attains its minimum. At small Vq [see the curve for Vq = 
0.01 in Fig. [IJa)], the values of A^ on the VK-stable branches approach the limit value, 
A^Townes ~ 5.85, which corresponds to the Townes soliton in the free 2D space 
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FIG. 1: (Color online) Chemical potential fi of the ground-state mode ("soliton") versus its norm 
N, at two fixed wavenumbers of the Penrose-tiling potential, A; = 1 (a) and k = 1.5 (b), and 
different values of its depth, Vq. Labels Cj and Aj {j = 1, 2, 3, 4) indicate VK-stable and unstable 
branches, respectively, while points Bj mark junctions between them. 

As said above, the main point in this work is the study of the solitons close to norm- 
minimizing points Bj, which are of obvious interest to the potential experiment. In Fig. [1] 
we observe that stable solitons with the minimum norm (i.e., smallest number of atoms) are, 
naturally, generated in the deepest potential, represented by families Al — Bl — CI. The 
norm attains its minimum, Nj^^^ = 1.304 (with /i = —0.518) at A; = 1 and Vq = 1 [point Bl in 
Fig. [D^a)]. We also observe that the stability range (the distance from the lower termination 
point to point Bl) in Fig. W^a) ior k = 1 is AA^ = A^(Cl) - A^(Bl) = 2.531 - 1.304 = 1. 227, 
which is ~ 20 times larger than AA^ = 5.510 — 5.460 = 0.05 in Fig. [T](b) for k = 1.5, at the 
same OL depth, Vq = 1. Generally, the comparison of Fig. [U^a) and Fig. [U^b) demonstrates 
that, for given Vq, the norm of the ground states strongly depends on the OL wavenumber, 
k. 
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FIG. 2: (Color online) The same as in Fig. [H but for the fixed depth of the OL potential, and 
different values of its wavenumber. 

The branch fi{N) with k = 0.9 in Fig. [2] is notably different from other branches with 
k > 1. Although a continuous dependence ^{N) is found in the range of C2 — B2b, no 
solutions have been found (the imaginary-time algorithm does not converge to them) between 
points B2b and B2a (the dashed segment B2b — B2a is depicted in Fig. |2]only as a guide to 
the eye). The algorithm again converges to the ground-state modes in the range of B2a — A2. 

Furthermore, a tail of segment B2b — C2 of this branch penetrates into the overcritical 
region, = 6.046 > A'xownes = 5.85. This feature is explained by the fact that the solitons 
found at < 0.9 (in particular, the ones marked by (3s, (3^, (3q in Fig. [2]) are actually double- 
humped structures, featuring pairs of spatially separated or almost fused density peaks [see 
Figs, m^c) and[5]^a), respectively]. 

Stability of the solitons. The VK criterion does not guarantee the full stability of solitons, 
as it does not capture instabilities associated with complex eigenvalues. To test the full 
stability, we simulated perturbed evolution of the solitons over a sufficiently long interval, 
typically t = 1000 (which covers, roughly, 10 diffraction times of the corresponding localized 
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states), adding small random perturbation to the initial conditions, with a relative amplitude 
~ 0.01. The modes whose evolution was tested in this way are indicated by arrows in Fig. |2l 
attached to symbols ai, and /Ss, /34, /Jg, which pertain to branches with k = \ and k = 0.9, 
respectively, and 75, that pertains to k = 1.1. The results of the evolution simulations are 
shown in Figs. [HE 

Figure |3] presents details of the stability test for the ground state on branch CI, marked 
by tti in Fig. [2] (for Vq = 1 and k = 1), with the norm and chemical potential = 2.098 
and /i = —1.027. This mode is stable. 

Figures m (a) and (b) display the evolution of the solitons taken near the junction points 
between the VK-stable and unstable segments of the fi{N) curves, for k = 1 and k = 0.9. 
Figure 111(a) pertains to the mode labeled 02 (with = 1) in Fig. [21 which evolves into the 
perturbed state depicted at t = 200 in Fig. 111(b). This mode is unstable, splitting into a 
set of density peaks located at different potential minima, which, however, do not tend to 
decay into dispersive waves. The evolution of another unstable mode, labeled by Pq in Fig. 
El is displayed in Figs. Hl^c) andlH^d). It splits into two parts, and eventually decays into 
the dispersive radiation either. 

Finally, Fig. [51 represents the perturbed evolution of the mode with larger norms (A^ > 4), 
for k = 0.9 and k = 1.1, which correspond to points and 75, respectively, labeled in Fig. 
121 In panels |5t^a,b) we again observe that the former (double-peak) mode, corresponding to 
k = 0.9, does not produce a stable soliton in the course of the perturbed evolution. However, 
Fig. El^d) demonstrates that the soliton corresponding to point 75 is stable. The eventual 
conclusion following from the analysis of the numerical results is that all the double-peak 
structures are unstable against splitting, irrespective of their formal compliance with the 
VK criterion, while the single-peak solitons are stable or not in the exact accordance with 
VK. 

Conclusion. We have studied the dynamics of 2D matter-wave solitons near the junction 
points between the stable and unstable branches of curves fi{N) for the soliton families 
supported by the interplay of the self-attractive nonlinearity and Penrose-tiling OL potential. 
These points are interesting to physical applications, as they correspond to the smallest 
number of atoms which is necessary to build 2D matter-wave solitons, or the smallest total 
power necessary for the making of spatial optical solitons. It was found that the shape and 
stability of such solitons crucially depend on the depth and period of the OL. A challenging 
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FIG. 3: (Color online.) The spatial structure of the stable localized mode supported by the 
quasiperiodic potential, labeled by ai in Fig. [2l (a) The Gaussian initial configuration ([3]) for 
Vb = 1 and k = 1, transformed by the imaginary-time relaxation into the ground state, which is 
shown in panel (b). Panel (c): The result of the perturbed evolution (in real time) at t = 1000. 
(d) The contour-plot profile of the underlying quasiperiodic potential with Vq = 1 and k = 1.0. 

problem is to extend the analysis to vortex solitons supported by quasi-periodic potentials jsj. 
This work was supported, in a part, by CONACyT/SEP 2012 and PROMEP CA Redes 
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FIG. 4: (Color online) The evolution of the stationary modes labeled by points 02 and /Sg in Fig. 
[21 for k = 0.9 and Vq = 1. (a) The shape of mode 02, with N = 1.304, /i = —0.518; (b) the result of 
the evolution at t = 1000. The final state is not a bound one, but it does not decay into radiation. 
(c,d) Mode (3^ with N = 1.319, ^ = —0.434, which splits into two parts, and eventually decays. 

projects (Mexico). 
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